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Complex dynamics of knotted filaments in shear flow 
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Abstract. - Coarse-grained simulations are used to demonstrate that knotted filaments in shear 
flow at zero Reynolds number exhibit remarkably rich dynamic behaviour. For stiff filaments that 
are weakly deformed by the shear forces, the knotted filaments rotate like rigid objects in the flow. 
But away from this regime the interplay between between shear forces and the flexibility of the 
filament leads to intricate regular and chaotic modes of motion that can be divided into distinct 
families. The set of accessible mode families depends to first order on a dimensionless number 
that relates the filament length, the elastic modulus, the friction per unit length and the shear 
rate. 



The interaction between a shear gradient and suspended 
objects can generate fascinating dynamical behaviour. G. 
Jeffrey [1] showed in 1922 that rigid bodies trace out com- 
plex periodic orbits that depend in detail on their shapes. 
For deformable objects [2] even richer behaviour is pos- 
sible. Red blood cells, for example, change shape with 
increasing shear rate [3]. Furthermore, the coupling to 
shear can lead to either unstable tumbling or to a steady 
state mode where the cells remain at a fixed angle to the 
flow while their outer membrane rotates like the tread- 
ing of a tank [4]. A similar crossover to tank treading is 
predicted for star polymers [5]. Experimental advances 
in single molecule techniques have made it possible to di- 
rectly observe the stretching and tumbling behaviour of 
individual DNA molecules [6]. This work inspired a great 
deal of theoretical research on the way that the polymer 
flexibility, shear and Brownian noise interact [7]. Sev- 
eral decades earlier it had been shown that filaments in 
the non-Brownian regime exhibit at least five distinguish- 
able regimes of motion as stiffness, length, and shear rate 
are varied [8]. These results are still the subject of ac- 
tive investigation by theorists [9]. The crossover between 
the Brownian and non-Brownian regimes has also been 
recently considered [10]. 

In this paper we use coarse-grained computer simula- 
tions to study the behaviour of knotted non-Brownian fil- 
aments in shear flow. Knots are a generic possibility for 
any long elastic objects and occur naturally in biologi- 



cally active DNA [11, 12]. In the limit of strong bending 
modulus A or weak shear rate 7 the knotted filament will 
take on its equilibrium shape [13], and rotate in a manner 
similar to that first predicted by Jeffrey [1]. Chiral knots 
should also migrate in the vorticity direction [14, 15], a 
hydrodynamic effect that has recently been observed for 
other objects including helical bacteria [16]. 

The focus of this paper, however, is what happens for 
stronger shear forces and/or for more flexible filaments, 
i.e. the regime where the knots can be tightened by the 
flow. To estimate where this crossover occurs we consider 
the following argument: in the limit of small local bond 
deformations, the bending energy of a knot of length L 
scales as ~ A/L so there is a force opposing knot tight- 
ening rsj A/L 2 . Neglecting the logarithmic factor [18], the 
drag on a slender filament in Stokes flow ~ nL. To tighten 
the knot, strands which are close to each other must be 
moved in opposite directions The typical velocity differ- 
ence will be ~ cry, where a is the filament width. Com- 
bining this with the drag and comparing to the bending 
force results in a dimensionless knot deformation number, 
a = A/ncrjL 3 . (At the crossover, the knot length « total 
filament length, so in calculating a we take L to be the 
total filament length.) This resembles the sperm number 
used with microscopic swimmers [19] but instead of deter- 
mining when the filament as a whole may be deformed, it 
indicates when the knot will be tightened. For larger a 
we expect the stiff knot regime, but for lower a the shear 
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Fig. 1: Predicted modes of behaviour as a function of the knot deformation number a. x is the flow direction, y is the direction of 
the shear gradient, and z is the vorticity direction. All configurations are for filaments with a 3i(— ) trefoil knot (a) Time series 
for filament in a II mode at a — 0. Time increases from left to right. The knot rotates clockwise around the z axis. The series 
spans approximately one period, the time between successive configurations is 1028. 7£o and migration is in the — z-direction. 
(b) Example configurations for different mode families over a range of a. Note that the orientation is different for the VII mode. 
Full animations of the modes are available online [17] (see Table 1 in the appendix). 



should cause significant deformation. 

Indeed, as illustrated in Fig. 1 for the case of the sim- 
ple trefoil knot, lowering a leads to a crossover from knots 
that remain close to their equilibrium shape, to a regime of 
surprisingly rich dynamical behaviour where the whole fil- 
ament exhibits intricate shape oscillations in time. These 
orbits may be grouped into a few distinct families com- 
prising very similar types of motion (modes). Some modes 
show regular, and others chaotic, motion. Different modes 
show distinct rates and directions of drift along the vor- 
ticity axis. A few families are accessible at each a, and 
the knot typically falls into one type of motion depending 
on initial conditions. Changing a therefore changes the 
modes that are accessible, as well as the probability that 
a certain mode is selected. In the appendix we list de- 
tails of videos, accessible online [17], which illustrate the 
modes. 

In the rest of this paper we describe how we simulate 
the knotted filaments and analyse the ensuing results. We 
apply a coarse-grained bead-spring model [20] . The inter- 
action potential between beads is [21]: 
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where ?i is the i tn bead position and f^ = r\ — fj. The 
first term in Eq. (1) allows the flexibility of the filaments to 
be varied by changing ft, the bending energy. The second 
term is excluded volume: H is the Heaviside step function 
which truncates the Lennard- Jones potential to be purely 



repulsive, a and e are the length and energy scales respec- 
tively. The third term is a FENE spring potential. We 
choose k = 30e/a 2 and Ro = 1.5a. 

We update bead positions using the Euler method [15]: 



fi(t + At) =f-(t) + 



v(ri) + ^2nij • fj 



At (2) 



where v(fi) is the applied flow, v = ^yx and fi is the 
force on the i th bead which results from the potential (1). 
Assuming zero Reynolds number, the hydrodynamic inter- 
actions, 1-Lij, can be approximated by the Rotne-Prager- 
Yamakawa interaction tensor [22] with viscosity 77, and 
taking the hydrodynamic radius as a/2. The natural input 
units for the simulation are <r, rj and /c, from which a natu- 
ral time unit of to = crrj/k follows. We used 7 = (150to) _1 , 
and verified that the stretching of individual bonds relaxes 
much more quickly than the characteristic time for shear 
induced motion and so should not influence the dynam- 
ics [23]. 

We simulated single filament rings with one knot, us- 
ing chains N = 50 beads (unless otherwise stated). We 
mainly study the chiral trefoil knot, standardly denoted 
3i [24], choosing, unless otherwise stated, the left-handed 
enantiomer (3i(— )), not the right-handed (3i(+)). 

Filaments were given a knotted configuration and then 
equilibrated for 6 x 10 5 to at finite temperature with no 
shear to generate random starting points. A 6 x 10 6 to ini- 
tialisation period was allowed before data was recorded for 
1.5 x 10 6 to- The relatively long initialisation was chosen to 
avoid transients. We ran 50 simulations for most parame- 
ter sets. The majority of runs exhibited the same mode of 
motion throughout the observation time, but for a small 
number (less than 1%) of runs a slightly longer initialisa- 
tion was necessary. The long transients may caused by the 
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system being close to a parameter value at which a mode 
appears/disappears [25]. Eq. (2) was typically integrated 
using At = 10 -2 although for some parameter choices, for 
example high ft, it was necessary to reduce this for nu- 
merical stability. We checked that using At = 10 ~ 3 gave 
equivalent results. 

Mode families were identified by visual inspection of 
their motion (see e.g. the videos in [17]) and by measuring 
their drift velocity in the vorticity direction, which distin- 
guished well between different modes. The identity was 
confirmed by two additional order parameters which mea- 
sure the direction of maximum extension and the asym- 
metry under a ir rotation about z. Detailed definitions, 
as well as plots of average values for individual runs, are 
given in the appendix. 

The elastic modulus of our bead-spring model A = 
kg [19], so that the knot deformation number defined ear- 
lier takes the form a = k/tjjL 3 . In Fig. 1 we consider 10 
different values of n corresponding to a = — 2.56 x 10 -3 . 
For the largest a, the filament remains in braid-like con- 
figurations [13] that characterise family VII in Fig. 1(b). 
The motion is composed partly of rotation of the config- 
uration and partly of tank-treading - a particular point 
moves around the contour. As flexibility is increased (a is 
lowered) there is a change from modes which rotate in the 
x-y-pldLiie to modes which have relatively large extensions 
in the z-direction. Interestingly, a similar shift was seen 
in experiment with linear filaments [8] . The first family in 
which the knot is significantly tightened is V. 

Some mode families show both regular and chaotic 
modes, sometimes at the same a, for example family II. 
Others showed only regular (VI), or only chaotic (IV), 
motion. To distinguish regular and chaotic modes, we cal- 
culated the largest Lyapunov exponent, g\ [26]: A second 
system was created with the bead positions rj each ran- 
domly displaced to Si and constrained so that d 2 = J2 i | 
Ti — Si | 2 /a 2 = 1. Both systems were integrated forward in 
time. After each 1.5 to, n — Si were rescaled by changing Si 
so as to make d — 1 . The Lyapunov exponent is then given 
by [26]: g\ = lim n _>oo0"n = lim n ^oo(2/3nt ) Yh=i m (4?')> 
where dj is the distance after the j th evolution. If the mea- 
sured <j n tended to zero or a positive constant as a func- 
tion of n the motion was identified as regular or chaotic 
respectively. This behaviour is illustrated in the inset of 
Fig. 2. 

Fig. 2 compares the movement of the average bead po- 
sition around its average drift in the z-direction for two 
modes from family II at a = 0, one regular and one 
chaotic. The regular mode simply oscillates. By con- 
trast the chaotic mode, whilst showing oscillations, also 
displays larger movements. The power-spectrum of the 
curve for chaotic motion, shown in the appendix, exhibits 
a power-law decay with an exponent of about minus 2 ( 
1.94 ±0.03), suggesting a random walk around the average 
drift. 

Average migration velocities were calculated by a lin- 
ear fits to the ^-displacements of the centre of resistance. 




x10 



Fig. 2: Comparing regular and chaotic runs with modes in 
family II at a = 0. The average ^-position of the filament 
beads with the average drift subtracted is plotted as a function 
of time. The inset shows the estimate of the largest Lyapunov 
exponent as function of time for the two runs, plotted on a 
log- log scale. 



Fig. 3 plots the averages over all runs, grouped into fam- 
ilies and then subdivided into regular and chaotic modes 
at each a. Error bars indicate the spread of velocities ob- 
served. For most they are smaller than the data points. 
Fig. 3 also shows that at some a there exist modes that 
migrate in opposite directions. The period of rotation of 
regular modes varies with a, approximately in the range 
2 — 6 x 10Ho. The shortest periods were seen for the high- 
est a. For example at a = the average period of modes 
in the II family was 6220 ± 40£o and that of family VII at 
a = 2.56 x 10" 3 was 2040±90t . Fig. 3 also shows the per- 
centage chance that a run with random initial conditions 
ends up in the particular mode family. 

We also considered simulations for other knot types at 
a = 0. For 3i(+) the same modes are seen but the migra- 
tion, and the orientation, is as expected, in the opposite z- 
direction. We see similar behaviour - regular and chaotic 
modes with migration - for more complex knots such as 
4i and 5i. For achiral 4i the distribution of migration ve- 
locities is symmetric about zero: all migrating modes have 
a partner with opposite migration direction. Depending 
on initial conditions, the 4i knot may thus migrate in the 
indirection, but that the average migration velocity over 
many runs would be zero. 

Finally, we consider the sensitivity of our results to 
changes in parameters and changes in simulation details. 
It should be kept in mind that for these dynamical sys- 
tems with behaviour that may depend sensitively on initial 
conditions, one would expect quantitative changes when 
simulation details are changed. The main thrust of our 
paper is qualitative, and so the most important tests will 
be whether the overall behaviour, i.e. the mode families, 
are robust to these changes. 

Firstly we consider the effect of hydrodynamic interac- 
tions by setting Hij = for i ^ j in Eq. (2). We find, 
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Fig. 3: Migration velocity in the ^-direction, averaged over all 
runs belonging to each family. The labels indicate families and 
are subdivided into regular and chaotic modes, denoted by the 
subscripts r and c respectively. Error bars show the standard 
deviations. The figures in brackets indicate the percentage of 
runs at each a that were observed to fall into each group. 



as expected, that the motion in the vorticity direction is 
a consequence of off-diagonal hydrodynamic interactions. 
Modes that resemble those of families V and VII were seen 
but none with large z-extensions. 

Secondly, linear filaments or unknotted rings with point 
force hydrodynamics simply align in the x-z plane with- 
out access to different shear velocities. Simulations must 
therefore explicitly take the finite thickness into account 
by considering the torque on individual beads [9]. We 
tested this sensitivity by using algorithms that include the 
torque, and find that, in contrast to unknotted filaments, 
similar mode families are observed. The knot forces the 
system out of the plane so that it always accesses different 
shear velocities, and this dominates. 

Thirdly, we checked how the behaviour is affected by 
changing ft, r] and 7 in such a way as to keep a fixed: in the 
absence discretisation effects, such changes of parameters 
should lead to descriptions of the same physical system 
and so the same behaviour is expected. We ran two sets 
of simulations where 7 was reduced by a factor of 10 and 
either n was decreased or 77 increased to compensate. We 
obtained very similar results. 

Finally, most of our results are for a fixed length fila- 
ment with N = 50 beads. It is interesting to investigate 
how sensitive our results are to the length L = Ncr. For 
example, for N = 50 significant tightening first occurs in 
family V at a = 0.64 x 10 -3 . We ran additional simula- 
tions with N = 40, 70 and 100 at a = 0.64 x 10~ 3 . In 
each case the majority of runs show a mode very similar 
to those in family V, see Fig. 4. The migration velocities 
are similar but decrease with L. The measured velocities 
are 1.60, 1.51, 1.35 and 1.27 xlO~ 4 a/to respectively (in 
each case all runs in the mode had exactly the same ve- 
locity to the accuracy given.) At other as we checked the 
results at different N were also qualitatively similar, al- 
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Fig. 4: Example configurations of the modes at a = 0.64 x 10 -3 
for N — 40, 50, 70 and 100. Full animations of the modes are 
available online [17] (see Table 1 in the appendix). 



though the agreement worsens at lower a. For example at 
a = 1.28 x 10 -3 modes like families VI and VII were seen 
for all N but for a < 0.16 x 10 ~ 3 we observed families 
for N = 100 that were qualitatively different to any seen 
for other N. In fact, we expect substantial differences at 
small a because the shape of the tighter knot is then fixed 
by the excluded volume of the chain and not just by the 
physics that enters into the derivation of a. In that regime, 
for fixed filament thickness, we expect the influence of the 
knot to become progressively smaller as L/cr — » 00. While 
it would be interesting to explore these effects further, at 
fixed 7, changing N while fixing a means that n must be 
increased as N 3 and the integration timestep correspond- 
ingly decreased for stability. Combining this with N 2 time 
for calculation of Hij gives a prohibitive ~ TV 5 scaling of 
simulation time. 

To summarise, we have demonstrated that knotted fil- 
aments in shear exhibit a rich dynamical behaviour with 
modes which can be divided into families. Some families 
have both regular and chaotic modes. Mode families mi- 
grate in different directions along the vorticity axis. The 
crossover from a stiff knot to the regime where multiple 
modes are possible can be described by a dimensionless 
number. In future work it may be interesting to consider 
more sophisticated treatments of the hydrodynamics that 
include effects such as lubrication. It may also be inter- 
esting to consider the effect of noise: Initial simulations 
suggest that fluctuations may alter the stability of modes 
leading to a variation or even flipping of migration velocity 
as function of noise strength. 

Experimentally, this behaviour would be most easily 
observable with macroscopic filaments in highly viscous 
solutions [8]. However, it may also be visible for DNA. 
For example, the P4 phage genome (common in knotting 
experiments [12]) is about 77 thermal persistence lengths 
long We estimate a crossover (a = 10 -3 ) at 7 « 2 x 10 3 s -1 
in water. The Weissenberg number « 10 so shear should 
be reasonably strong compared to thermal effects. 

Appendix. — We present a range of animations of 
modes from the families described in the main text. The 
filenames, along with additional information, are listed in 
Table 1. All animations are of duration 9000£o and of 
simulations at a shear rate of 7 = (150£q) -1 with a 3i(— ) 
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knot. The green sections of the filaments are markers to 
allow the motion to be more easily followed. 



Table 1: Examples of animations of modes belonging to the 
various families described in the main text. Files may be ac- 
cessed online [17]. Animations were created using VMD [27]. 
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We next briefly discuss the two order parameters that 
were used to help group runs into mode families. The first, 
(f). was the angle of the direction of maximum extension 
to the z-axis, allowed to vary between and ir/2. <p was 
determined by finding the eigenvector of the largest eigen- 
value of the radius of gyration tensor. The second, C2, 
was defined as follows 



C 2 



1 
~NR 



y^min(\fj 



(3) 



where R is the average bead separation and r- are the 
bead positions rotated about the z-axis by ir in the centre 
of mass frame: the minimum distance from each bead to 
a bead in the rotated configuration is summed. Smaller 
values of C2 indicate configurations which are closer to 
being symmetric under a tt rotation. 

Figs. 5 (a) and (b) show the values for these two order 
parameters for different a for the N = 50 results. Each 
point is the average over one of fifty runs - they are plotted 
in an arbitrary order. It should be emphasised that all the 
points within two consecutive vertical lines are for different 
runs for the same a - the different positions along the x- 
axis within each section are irrelevant. 

We also include a plot of the power-spectrum of the 
data for the chaotic mode plotted in Fig. 6. This was 
obtained by taking the modulus-squared of the discrete 
Fourier transform of the displacement of the average bead 
position around its overall drift. As may be seen from 
Fig. 6, the exponent of the decay is close to -2 (the mea- 
sured value is -1.94 ± 0.03). 



Fig. 5: The values of the order parameters, averaged over sin- 
gle runs for N — 50 filaments at different values of a. (a) The 
angle of the direction of maximum extension to the z-axis, 0. 
The averages of (f) for each run are plotted for a given a in an 
arbitrary order. The labels indicate the modes to which the 
different groups of points correspond, (b) The same as (a) but 
for C2, an order parameter to detect two- fold symmetry about 
the z-axis. Lower values indicate more symmetric configura- 
tions. It should be emphasised that all the points within two 
consecutive vertical lines correspond to different runs at the 
same a - the positions along the x-axis within each section are 
irrelevant. 




Fig. 6: Power spectrum calculated by discrete Fourier trans- 
form of the displacement about the average drift for the chaotic 
mode plotted in Fig. 2. The dashed line has a slope of -2. 
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